#include "ppdef.h"


program dineof
 use ufileformat
 use initfile
 implicit none

#ifdef PARALLEL_ARPACK
include "mpif.h"
#endif


 integer              ::   maxm, maxn, maxnev, maxncv, ldv, ldu,nev,ncv,&
      &                        M,N,I,J,IMISS,IMEAN,ICROSS,IMISST,IM,k,IT,P,NN,&
      &                        NITE,ipid, ierror,R,nitemax,q,getpid,&
      &                        MAKENIT,ival,rec,nbvars,eofno,testN,norma

 integer,pointer      :: imax(:),jmax(:),first(:)

! parameter               (maxnev=1000, maxncv=1010)

 real, allocatable    :: meanA(:),& 
      &                        w(:), V(:,:), U(:,:),&
      &                        workl(:), workd(:),&
      &                        S(:,:), resid(:),xcr(:,:),X2(:,:)
 real                 :: tol, t, sumSV, XMEAN, VALEX,Y,DIFF,VAL,valexmask
 real                 :: t1,t2,t3,mytol,testmean,sumVarIni,SumVarEnd
 real, pointer        :: X(:,:),mask(:,:),Xp(:,:),fileMean(:),fileStd(:)

 ! indexes of cross-validation points choosen by user
 real, pointer        :: cvp_indexes(:,:)
 real                 :: valexc

 real, allocatable    :: XEX(:,:),VALC(:),sst(:,:,:),img(:,:,:)
 integer, allocatable :: IEX(:),JEX(:)

 logical,allocatable  :: select(:)
 logical              :: istrans
 character*72            comment

 character(len=200),pointer   :: sstfile(:),maskfile(:),resultfnames(:)
 character(len=200)   :: DirOutput, initfilename, cloudsfname,               &
      logfile

 integer, pointer     :: seed(:)
 integer              :: iargc

 integer              :: istat

#ifdef DIAG_CROSSVAL
 character(len=200)   :: DirDiag
 character(len=200),                                                          &
      pointer         :: diagfnames(:)

#endif

#ifdef PARALLEL_ARPACK
integer              :: comm, ierr,myrank,nbprocs
#else
integer,parameter    :: myrank=0     
#endif

 include 'ReadMatrix.h'


 interface
   subroutine valsvd(U,S,V,IEX,JEX,VAL,M,N,maxnev,maxncv,ldu,ldv,NN)
    integer,intent(in) :: maxnev,maxncv,ldu,ldv
    real, intent(in) :: U(ldu, maxnev),V(ldv,maxncv),S(maxncv,2)
    integer,intent(in) :: M,N,NN,IEX,JEX
    real,intent(out) :: VAL
   end subroutine valsvd
 end interface


 interface
   subroutine writeMatrix(X,xmean,valex,sstfile,maskfile,norma,M,N,imax,jmax,first,DirOutput,valc,fileMean,fileStd)
    real,intent(in)                  :: X(:,:),fileMean(:),fileStd(:)
    integer,intent(in)               :: N,M,imax(:),jmax(:),first(:)
    real,intent(in)                  :: xmean,valc(:)
    real,intent(inout)               :: valex
    integer,intent(in)               :: norma
    character(len=200),intent(in)    :: sstfile(:),maskfile(:)
    character(len=200),intent(in)    :: DirOutput
   end subroutine writeMatrix
 end interface

 interface
   subroutine writeSVD(DirOutput,s,nev,u,n,m,v,sumVarIni,sumVarEnd)
    character(len=200),intent(in)     :: DirOutput
    real,intent(in)                   :: s(:,:),u(:,:),v(:,:),sumVarIni,sumVarEnd
    integer,intent(in)                :: nev,n,m
   end subroutine writeSVD
 end interface

 !      ----------------------------------------
 !      Initiate MPI sequence
 !      ----------------------------------------  

#ifdef PARALLEL_ARPACK
 call MPI_INIT( ierr )
 comm = MPI_COMM_WORLD
 CALL MPI_COMM_RANK( comm, myrank, ierr )
 CALL MPI_COMM_SIZE( comm, nbprocs, ierr )
 
#endif

 !      ----------------------------------------
 !      THIS DOES NOT WORK WITH MPI: check if DINEOF has been called properly
 !      ----------------------------------------  

! if(iargc().ne.1) then
!   write(stderr,*) 'DINEOF usage:'
!   write(stderr,*) 'dineof <init filename> ',iargc(),myrank,nbprocs,comm
!   stop
! end if

! call getarg(1,initfilename)

 !-------------------------------------------------------------------------
 !      Read numerical parameters from file initfilename
 !-------------------------------------------------------------------------
 initfilename='dineof.init'
 call getInitValue(initfilename,'nev',nev) 
 call getInitValue(initfilename,'ncv',ncv) 
 call getInitValue(initfilename,'tol',tol) 
 call getInitValue(initfilename,'rec',rec) 
 call getInitValue(initfilename,'eof',eofno) 
 call getInitValue(initfilename,'norm',norma) 
 call getInitValue(initfilename,'Output',DirOutput) 
#ifdef DIAG_CROSSVAL
 call getInitValue(initfilename,'DirDiag',DirDiag) 
 call getInitValue(initfilename,'diags',diagfnames) 
#endif
 call getInitValue(initfilename,'data',sstfile) 
 call getInitValue(initfilename,'mask',maskfile) 
 call getInitValue(initfilename,'results',resultfnames) 
 call getInitValue(initfilename,'seed',ipid) 

 if (stdout.ne.6) then
   call getInitValue(initfilename,'logfile',logfile,default='dineof.out') 
   open(stdout,file=logfile,status='unknown')
 end if

 write(stdout,*) '****************************************************'
 write(stdout,*) 'Numerical data read'
 write(stdout,*) 'You entered the values:'
 write(stdout,*) 'number of EOF modes you want to compute',nev
 write(stdout,*) 'maximal size for the Krylov subspace',ncv
 write(stdout,*)

 call flush(stdout,istat)

 if(norma.eq.1) then
   write(stdout,*) 'You asked for the normalisation of the input matrices'
 else
   write(stdout,*) 'You asked not to normalise the input matrices'
 end if

 write(stdout,*)

 if(eofno.eq.1) then
   write(stdout,*) 'The right and left EOFs will be written in directory ',trim(DirOutput)
 else
   write(stdout,*) 'You asked not to write the left and right EOFs used for the reconstruction'
   write(stdout,*) '(to set this option ON, change eofno=0 in file Input/numerical.data)'
 end if
 write(stdout,*)

 write(stdout,*) '****************************************************'
 !-----------------------------------

 maxnev=nev
 maxncv=ncv
 nitemax=35

 allocate(select(maxncv))

 allocate(VALC(0:nev+1))

 nbvars = size(sstfile)

 do q=1,nbvars
   write(stdout,*)'You entered filenames ', trim(sstfile(q)) 
   write(stdout,*)'                      ', trim(maskfile(q)) 
 end do


 !--------ReadMatrix-------------------------------------------------------------------
 !     Reads input matrices and the exclusion value VALEX 
 !     Reshapes the input matrices into one matrix with size M*N (spatial and temporal)
 !     Input(sstfile(s), maskfile(s),norma)
 !     Output(X,M,N,VALEX,mask,imax,jmax,first,fileMean,fileStd)
 !-------------------------------------------------------------------------------------

 call ReadMatrix(sstfile,X, M, N,VALEX,maskfile,norma,imax,jmax,first,fileMean,fileStd)

 !--------------------------------------------
 write(stdout,*) '***'
 write(stdout,*) 'Matrix loaded '
 write(stdout,*) 'Size of the matrix = ', m,' by ', n
 write(stdout,*) '***'
 !--------------------------------------------

 maxm = M
 maxn=N

 ldu = maxm
 ldv=maxn

 


 ! Transpose matrix if M le N----------------------------------------------
 istrans=.false. 
 if (M.lt.N) then

    allocate(X2(N,M))
    X2=transpose(X)
    deallocate(X)
    allocate(X(N,M))
    X=X2
    deallocate(X2)
  
    I=M 
    M=N 
    N=I 
    maxm = M
    maxn=N
    ldu = maxm
    ldv=maxn
    allocate(meanA(maxm),& 
      w(ldu), V(ldv,maxncv), U(ldu, maxnev),&
      workl(maxncv*(maxncv+8)), workd(3*maxn),&
      S(maxncv,2), resid(maxn))

    ISTRANS=.true. 
!   stop "PLEASE PROVIDE TRANSPOSED MATRIX SINCE M GT N" 
 else


    allocate(meanA(maxm),& 
      w(ldu), V(ldv,maxncv), U(ldu, maxnev),&
      workl(maxncv*(maxncv+8)), workd(3*maxn),&
      S(maxncv,2), resid(maxn))
 end if 


 !----------------------------------------------------------------
 !--------------------------------------------------------

 !     ----------------------------------------------------------
 !     First statistics on reshaped matrix
 !     Calculate mean, missing data, and cross validation points
 !     ---------------------------------------------------------- 

 XMEAN=0 
 IMEAN=0 

 do J=1,N 
   do I=1,M
     if(X(I,J).ne.VALEX) then 
       IMEAN=IMEAN+1 
       XMEAN=XMEAN+X(I,J) 
     end if
   end do
 end do

 XMEAN=XMEAN/FLOAT(IMEAN) 
 IMISS=M*N-IMEAN 

 write(stdout,*) 'Missing data:',IMISS,'out of ', M*N,' (',IMISS*100./(M*N),'%)'
 write(stdout,*)

 if (.not.presentInitValue(initfilename,'clouds')) then
   ICROSS=min(0.01*M*N+40,0.03*M*N)
 else
   ! user chooses the cross-valiation points by the entry "clouds"
   call getInitValue(initfilename,'clouds',cloudsfname)
   call uload(cloudsfname,cvp_indexes,valexc)
   write(stdout,*) 'User-specified cross-validation points: ',trim(cloudsfname)
   icross = size(cvp_indexes,1)
 end if

 write(stdout,*) 'Number of cross validation points ',ICROSS 
 write(stdout,*) '****************************************************************'
 write(stdout,*)

 call flush(stdout,istat)

 IMISST=IMISS+ICROSS

 if(IMISST.gt.(M*N)) stop'Not enough data' 

 !     ------------------------------------------------------
 !     Declare and initialise IEX(IMISST),JEX(IMISST)
 !     data actually mising from the input matrix
 !     ------------------------------------------------------ 

 allocate(XEX(IMISST,0:nev),IEX(IMISST),JEX(IMISST))


 ! M*N OPS 
 IM=1
 do J=1,N 
   do I=1,M 
     if(X(I,J).eq.VALEX) then 
       IEX(IM)=I 
       JEX(IM)=J            
       IM=IM+1 

     else 
       !     -----------------
       !     substract mean 
       !     -----------------
       X(I,J)=X(I,J)-XMEAN 
     end if
   end do
 end do

 if(IM-1.ne.IMISS) stop '???' 

 !initial variance of the matrix X
 sumVarIni = M * N * sum(X**2,X.ne.valex)/count(X.ne.valex)
 write(stdout,*)'number of points ',M*N
 write(stdout,*)'number of points non missing',count(X.ne.valex)
 write(stdout,*)'energy ',sum(X**2,X.ne.valex)
 write(stdout,*)'varini ',sum(X**2,X.ne.valex)/count(X.ne.valex)


 !-------------------------------------------------------------------------------
 ! Add pseudo missing data points for cross validation, but save the values into 
 ! an array 
 ! declare XEX(IMISST,0:N) 
 ! XEX(II,NN) contains the interpolation at point IEX(II),JEX(II) when NN modes are used 
 ! For NN=O initial data are stored 



 if (.not.presentInitValue(initfilename,'clouds')) then
   ! choose cross-validation points internally

   !----random data set (random cross-validation points) ----------------------------------

   ! ipid=getpid() 
   call random_seed(SIZE = R)
   allocate(seed(R))
   seed = ipid
   call random_seed(put = seed(1:R))
   deallocate(seed)

   K=0 

   do while (k.lt.icross)

     call random_number(Y)
     I=floor(Y*M)+1
     call random_number(Y)  
     J=floor(Y*N)+1 



     if(X(I,J).ne.VALEX) then 
       ! for cross validation

       K=K+1 
       IEX(IMISS+K)=I 
       JEX(IMISS+K)=J 
       XEX(IMISS+K,0)=X(I,J)
       X(i,j) = valex 
     end if

   end do

 else
   ! user chooses the cross-valiation points by the entry "clouds"
   do k=1,icross
     i = cvp_indexes(k,1)
     j = cvp_indexes(k,2)

     if(X(I,J).ne.VALEX) then 

       IEX(IMISS+K) = i
       JEX(IMISS+K) = j
       XEX(IMISS+K,0)=X(I,J)

       ! for debuging
       !if (k.le.10) write(stdout,*) 'first cross validation point ',X(i,j)+XMEAN

       X(i,j) = valex
     else
       write(stderr,*) 'Error in ',trim(cloudsfname),X(I,J)
       stop
     end if
   end do

   deallocate(cvp_indexes)
 end if

 !initialise missing values (=VALEX) to zero
 do J=1,N 
   do I=1,M
     if(X(I,J).eq.VALEX) then 
       X(I,J)=0
     end if
   end do
 end do

 !-------------------------------------------------------------
 !    Estimates the total variance

 call svariExp&
      &     (X, M, N, maxm, maxn, sumSV, w)
 !------------------------------------------------------------

 !-------------------------------------------------------------
 ! Now first guess of EOF 1 U,V,D can be calculated from XCR 
 !     -----------
 !-------------------------------------------------------------
 ! call LANCZOS 



 MAKENIT=1
 mytol=0

 call cpu_time(t1)
 IT=0
 call ssvd_lancz(&
#ifdef PARALLEL_ARPACK
      &       comm, &
#endif
      &       X, M, N, maxm, maxn, nev, ncv, maxnev, maxncv, ldv,& 
      &       ldu, w, V, U, workl, workd, S, resid, select, tol, IT,&
      &       MAKENIT,mytol)

 call cpu_time(t2)

 write(stdout,*)'time in lanczos for 1st cycle in DINEOF',t2-t1
 write(stdout,*)

 ! VALSVD: should give value at point I,J for P modes 
 ! Fill in missing data points with first guess

 do K=1,IMISST

   NN=1

   call valsvd(U,S,V,IEX(K),JEX(K),VAL,M,N,maxnev,maxncv,ldu,ldv,NN) 
   X(IEX(K),JEX(K))=VAL

 end do



 ! now loop on 'nev' modes asked for the reconstruction 
 !--------------------------------------------------------------------
 !--------------------------------------------------------------------
 ! 

 valc(0)=0

 do P=1,nev 

   ! Iterations.

   write(stdout,*) 'computing eigenvector',P,' (of',nev,')'
   call flush(stdout,istat)

   do NITE=1,nitemax


     !test no exploit of precedent iterations:it=0
     IT=1

     call ssvd_lancz(&
#ifdef PARALLEL_ARPACK
          &       comm, &
#endif
          &       X, M, N, maxm, maxn, P, ncv,maxnev, maxncv, ldv,& 
          &       ldu, w, V, U, workl, workd, S, resid, select, tol, IT,MAKENIT,mytol)


     do K=1,IMISST 
       NN=P 
       call valsvd(U,S,V,IEX(K),JEX(K),VAL,M,N,maxnev,maxncv,ldu,ldv,NN) 
       X(IEX(K),JEX(K))=VAL 
     end do

   end do


   ! Store best guess with P EOFs: 
   do K=1,IMISST
     XEX(K,P)=X(IEX(K),JEX(K)) 
   end do

   ! Calculate Cross validator 
   VAL=0 

   do K=1,ICROSS 
     DIFF=X(IEX(IMISS+K),JEX(IMISS+K))-XEX(IMISS+K,0) 
     VAL=VAL+DIFF*DIFF 
   end do

   VALC(P)=sqrt(VAL/FLOAT(ICROSS)) 

   mytol=(1./100.)*(valc(p)-valc(p-1))/valc(p)

 end do
 call cpu_time(t3)




 !-----------------------------------------------------------------
 ! Finished all values for P 
 !-----------------------------------------------------------------
 ! Find value of P which corresponds to the minimum value of VALC 

 P=1 
 VAL=VALC(1) 
 do I=2,nev 
   if(VALC(I).lt.VAL) then 
     VAL=VALC(I) 
     P=I 
   end if
 end do

#ifdef DIAG_CROSSVAL
 ! TO DO!!!
 !transpose is not implemented in this ifdef section!

 ! Recalculate the EOFs for optimal number of modes (P) without the cross-validation points 
 ! for diagnostic purposes

 allocate(Xp(M,N))

 do NITE=1,nitemax

   call ssvd_lancz(&
#ifdef PARALLEL_ARPACK
      &         comm, &
#endif
        &       X, M, N, maxm, maxn, P, ncv, maxnev, maxncv, ldv,& 
        &       ldu, w, V, U, workl, workd, S, resid, select, tol, IT,MAKENIT,mytol)

   do K=1,IMISS 
     NN=P 
     call valsvd(U,S,V,IEX(K),JEX(K),VAL,M,N,maxnev,maxncv,ldu,ldv,NN) 
     X(IEX(K),JEX(K))=VAL 
   end do


 end do

   ! ReCalculate Cross validator 
   VAL=0 

   do K=1,ICROSS 
     DIFF=X(IEX(IMISS+K),JEX(IMISS+K))-XEX(IMISS+K,0) 
     VAL=VAL+DIFF*DIFF 
   end do

   VAL=sqrt(VAL/FLOAT(ICROSS)) 

   write(stdout,*) 'check cross-validator ',val


 if(rec.eq.1) then
   do i=1,M
     do j=1,N
       Xp(i,j)=Xmean
       do k=1,P
         Xp(i,j)=Xp(i,j)+U(i,k)*S(k,1)*V(j,k)
       end do
     end do
   end do
 end if

 call writeSVD(DirDiag,s,nev,u,n,m,v,sumVarIni,SumVarEnd)

 call writeMatrix(Xp,xmean,valex,diagfnames,maskfile,norma,M,N,imax,jmax,first,DirDiag,valc,fileMean,fileStd)

#endif

 !--------------------------------------------------------------

 ! Now optimal number P of EOFs is known: 
 ! Retrieve best guess from that number 

 do K=1,IMISS 
   X(IEX(K),JEX(K))=XEX(K,P) 
 end do


 !--------------------------------------------------------------
 ! Use now the points put aside from initial data set 

 do K=1,ICROSS 
   X(IEX(IMISS+K),JEX(IMISS+K))=XEX(IMISS+K,0) 
 end do




 IT=1
 ! Make last iterations from there but only on IMISS, not IMISST 

 do NITE=1,nitemax

   call ssvd_lancz(&
#ifdef PARALLEL_ARPACK
      &         comm, &
#endif
        &       X, M, N, maxm, maxn, P, ncv, maxnev, maxncv, ldv,& 
        &       ldu, w, V, U, workl, workd, S, resid, select, tol, IT,MAKENIT,mytol)

   do K=1,IMISS 
     NN=P 
     call valsvd(U,S,V,IEX(K),JEX(K),VAL,M,N,maxnev,maxncv,ldu,ldv,NN) 
     X(IEX(K),JEX(K))=VAL 
   end do

 end do

 SumVarEnd =  M * N * sum(X**2,X.ne.valex)/count(X.ne.valex)

 ! Reconstruct the matrix with the EOFs retained as optimal
 if (rec.eq.1) then
   do i=1,M
     do j=1,N
       X(i,j)=0
       do k=1,P
         X(i,j)=X(i,j)+U(i,k)*S(k,1)*V(j,k)
       end do
     end do
   end do
 end if

 ! Almost finished: Retranspose if necessary 
 ! If(ITRANS) THEN 
 ! CALL ??TRANS??(X,M,N) 
 ! I=M 
 ! M=N 
 ! N=I 
 ! ENDIF 


 write(stdout,*)
 write(stdout,*)'DINEOF finished!'
 write(stdout,*)
 write(stdout,*)'number of eigenvalues retained for the reconstruction:',P
 write(stdout,*)'expected error calculated by cross-validation',valc(P)
 write(stdout,*)'total time in lanczos process', t3-t1,'seconds'
 write(stdout,*)
 write(stdout,*)'Now writing data...'

if(myrank.eq.0) then
  if(eofno.eq.1) then
      if (istrans) then

          ! u and v have to be swaped if matrix has been transposed
          call writeSVD(DirOutput,s,P,v,m,n,u,sumVarIni,SumVarEnd)
      else
          call writeSVD(DirOutput,s,P,u,n,m,v,sumVarIni,SumVarEnd)
     endif
  endif
endif

if (istrans) then

 allocate(X2(N,M))
 X2=transpose(X)
 deallocate(X)
 allocate(X(N,M))
 X=X2
 deallocate(X2)
    I=M 
    M=N 
    N=I 
    maxm = M
    maxn=N
endif

 where (X.ne.valex) X=X+xmean 

if(myrank.eq.0) then
 call writeMatrix(X,xmean,valex,resultfnames,maskfile,norma,M,N,imax,jmax,first,DirOutput,valc,fileMean,fileStd)
endif
 
 deallocate(X,imax,jmax,first,fileMean,fileStd,sstfile,maskfile,resultfnames)

 write(stdout,*)
 write(stdout,*)'...done!'
 write(stdout,*)

 if (stdout.ne.6) close(stdout)

#ifdef PARALLEL_ARPACK
call MPI_FINALIZE(ierr)
#endif

end program dineof






